Results for Supplemental Material for Hitchcock & Frank 2023
Peter Frank Hitchcock - LNCC Lab, Brown University
2024-04-05
- Load data from studies 1 and 2
- Plotting delay
- Supplemental
material model validation (although learning curves exp 2 plot is in
key_results) - Parameter recovery
- Simulation with a high difference in phi
- Supplementary Figure 7
- PST/Stimulus valuation phase analyses
- 3/19/24 - ethnicity as per Ethnicity.identified var from Prolific from shared drive (shared drive doc has other demographic tallies)
- Questionnaires
sapply(c(
"rjson",
"data.table",
"dplyr",
"ggplot2",
"stringr",
"purrr",
"foreach",
"doParallel",
"patchwork",
"lme4",
"lmerTest",
"testit",
"latex2exp",
"ggpubr"
),
require, character=TRUE)## Loading required package: rjson## Loading required package: data.table## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: ggplot2## Loading required package: stringr## Loading required package: purrr##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: doParallel## Loading required package: iterators## Loading required package: parallel## Loading required package: patchwork## Loading required package: lme4## Loading required package: Matrix## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Loading required package: testit## Loading required package: latex2exp## Loading required package: ggpubr## rjson data.table dplyr ggplot2 stringr purrr foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## doParallel patchwork lme4 lmerTest testit latex2exp ggpubr
## TRUE TRUE TRUE TRUE TRUE TRUE TRUEsf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()## ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value ? ?
## visible FALSE FALSE
## ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value ? ?
## visible FALSE FALSE
## ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearner.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/aModelUtils.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModel.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/RecodeDfs.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayQVals.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value ?
## visible FALSEDefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))Load data from studies 1 and 2
s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv")
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv")
s1_pst <- read.csv("../data/cleaned_files/s1_PST_deident.csv") %>% rename(ID=deident_ID)
s2_pst <- read.csv("../data/cleaned_files/s2_PST_deident.csv") %>% rename(ID=deident_ID)Harmonize variables and create some separate vars
# Study 2 harmonize
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive"
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2
s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))
s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))Define paths and functions
# MLE results
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function
rm <- function(path, model_str) read.csv(paste0(path, model_str))Plotting delay
s1_train_delay_summs <- s1_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.s1_tr_delay_summs <-
Rmisc::summarySEwithin(s1_train_delay_summs,
measurevar = "m",
withinvars = c("condition", "delay_trunc"),
idvar = "ID")## Automatically converting the following non-factors to factors: condition, delay_truncs2_train_delay_summs <- s2_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.s2_tr_delay_summs <-
Rmisc::summarySEwithin(s2_train_delay_summs,
measurevar = "m",
withinvars = c("condition", "delay_trunc"),
idvar = "ID")## Automatically converting the following non-factors to factors: condition, delay_trunca <- ggplot(s1_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) +
geom_line(size=2, position = position_dodge(width = .2)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") +
geom_hline(yintercept = .5, size=2, color="black") +
geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
ga + ap + ft + tol + ylim(.48, .82) + tp +
scale_color_manual(values=c("purple", "orange"),
labels=c("cognitive", "overt")) + ylab("Proportion correct") + xlab("Delay") +
ggtitle("Experiment 1")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.b <- ggplot(s2_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) +
geom_line(size=2, position = position_dodge(width = .2)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") +
geom_hline(yintercept = .5, size=2, color="black") +
geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
ga + ap + ft + lp + ylim(.48, .82) + tp +
scale_color_manual(values=c("purple", "orange"),
labels=c("cognitive", "overt")) + ylab("") + xlab("Delay") +
ggtitle("Experiment 2")delay_comb_plot <- a+bdelay_comb_plot
#ggsave("../paper/figs/fig_parts/delay_plot.png", delay_comb_plot, height=7, width=10, dpi=700)Supplemental material model validation (although learning curves exp
2 plot is in key_results)
6/5/23 — reran sims to mkae sure completely updated/bug fixed version of par results used
s1_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")
s1_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")Capture of delay effects
out_a <- SimplePlotSimTrainingDelay(emp_df=s1_train, sim_df=s1_train_sim_m1_eb)## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## Warning: Duplicated aesthetics after name standardisation: colourout_b <- SimplePlotSimTrainingDelay(emp_df=s2_train, sim_df=s2_train_sim_m1_eb) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## Warning: Duplicated aesthetics after name standardisation: colourout_a## Warning: Duplicated aesthetics after name standardisation: colour
out_b## Warning: Duplicated aesthetics after name standardisation: colour
# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_1.png", out_a, height=8, width=12, dpi=700)
# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_2.png", out_b, height=8, width=12, dpi=700)Worst-option Proportions in Test Phase and comparison to same model without CK
worst_plots_a_ck <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")w_a_ck <- worst_plots_a_ck[[1]]+worst_plots_a_ck[[2]]w_a_ck 
worst_plots_b_ck <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")w_b_ck <- worst_plots_b_ck[[1]]+worst_plots_b_ck[[2]]w_b_ck## Warning: Removed 1 rows containing missing values (`geom_point()`).
# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_1.png", w_a_ck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_2.png", w_b_ck, height=6, width=15, dpi=700)Comparative worst plots
s1_test_sim_m1_eb_no_ck <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK16278.csv")
s2_test_sim_m1_eb_no_ck <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK31508.csv")worst_plots_nck_a <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")w_a_nck <- worst_plots_nck_a[[1]]+worst_plots_nck_a[[2]]w_a_nck
worst_plots_nck_b <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")Captures the mis-spec
w_b_nck <- worst_plots_nck_b[[1]]+worst_plots_nck_b[[2]]w_b_nck
# ggsave("../paper/figs/fig_parts/NOCK_simvsemp_worst_plot_1.png", w_a_nck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/NOCKsimvsemp_worst_plot_2.png", w_b_nck, height=6, width=15, dpi=700)Parameter recovery
m1_par_recov <-
read.csv("../model_res/par_recov_clean/pr_opts_mle/best_pr/par_recov_BEST_EMPIRICAL_BAYES_mv_gaussstudy_1_train_SIT_RunMQLearnerDiffDecayToPessPrior_76217.csv") # Confirmed created 3/28 after the 3/23 bug fix on decay eb_recov_pars_m1 <- m1_par_recov[grep("EB_recovered", names(m1_par_recov))]
simmed_pars_m1 <- m1_par_recov[grep("simmed", names(m1_par_recov))]Phi difference
this_eb_pr_df <-
data.frame("simulated"=simmed_pars_m1$cog_phi_simmed-simmed_pars_m1$overt_phi_simmed,
"EB_recovered"=eb_recov_pars_m1$cog_phi_EB_recovered-eb_recov_pars_m1$overt_phi_EB_recovered)
phi_diff <- ggplot(this_eb_pr_df, aes(x=simulated, y=EB_recovered)) +
geom_line(aes(x=simulated, y=simulated), size=4) +
geom_point(size=8, pch=21, fill="gray57", alpha=.8) +
stat_cor(method="spearman", size=6, label.y=.9) +
ggtitle(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
ga + ap + tp +
ylab("Recovered") + xlab("Simulated")
phi_diff
#ggsave("../paper/figs/fig_parts/pr/pr_phi_diff.png", phi_diff, height=4, width=8, dpi=700)Percent correctly recovering the sign difference
this_eb_pr_df$sim_sign <- if_else(this_eb_pr_df$simulated >= 0, 1, 0)
this_eb_pr_df$eb_sign <- if_else(this_eb_pr_df$EB_recovered >= 0, 1, 0)
# this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign
# (this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1
ta <- table((this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1)
ta[2]/sum(ta)## 1
## 0.7142857plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$epsilon_simmed, eb_recov_pars_m1$epsilon_EB_recovered,
"$\\epsilon", stat_pos = .32)
plot
# ggsave("../paper/figs/fig_parts/pr/eps.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$cog_phi_simmed, eb_recov_pars_m1$cog_phi_EB_recovered,
"$\\phi^{Cog}")
plot
# ggsave("../paper/figs/fig_parts/pr/phi_cog.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$overt_phi_simmed, eb_recov_pars_m1$overt_phi_EB_recovered,
"$\\phi^{Overt}")
plot
# ggsave("../paper/figs/fig_parts/pr/phi_overt.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$choice_LR_simmed, eb_recov_pars_m1$choice_LR_EB_recovered,
"$\\alpha_{CK}")
plot
# ggsave("../paper/figs/fig_parts/pr/ck.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$explor_scalar_simmed, eb_recov_pars_m1$explor_scalar_EB_recovered,
"ES")
plot
# ggsave("../paper/figs/fig_parts/pr/es.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$q_LR_simmed, eb_recov_pars_m1$q_LR_EB_recovered,
"$\\alpha_{Q}")
plot
# ggsave("../paper/figs/fig_parts/pr/qlr.png",
# plot, height=4, width=8, dpi=700)plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$beta_simmed, eb_recov_pars_m1$beta_EB_recovered,
"$\\beta", stat_pos = 25)
plot
Simulation with a high difference in phi
Load model fits
m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")
m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")
m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Find a participant with relatively high diff in decay
m1_study1_eb$cog_phi[25]## [1] 0.6817201m1_study1_eb$overt_phi[25]## [1] 0.1300788m1_study1_eb[25, ]## # A tibble: 1 × 12
## # Groups: ID [1]
## X epsilon q_LR cog_phi overt_phi choice_LR explor_scalar beta
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 25 0.0365 0.732 0.682 0.130 0.0425 0.220 4.97
## # ℹ 4 more variables: convergence <int>, nll <dbl>, AIC <dbl>, ID <int>This did correspond to a big difference in performance empirically
s1_train %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.## # A tibble: 8 × 3
## # Groups: valence_and_probability [4]
## valence_and_probability condition m
## <chr> <chr> <dbl>
## 1 40-10_punishment cognitive 0.6
## 2 40-10_punishment overt 0.95
## 3 40-10_reward cognitive 0.9
## 4 40-10_reward overt 1
## 5 90-10_punishment cognitive 0.525
## 6 90-10_punishment overt 0.775
## 7 90-10_reward cognitive 0.05
## 8 90-10_reward overt 1s1_sit %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.## # A tibble: 8 × 3
## # Groups: valence_and_probability [4]
## valence_and_probability condition m
## <chr> <chr> <dbl>
## 1 40-10_punishment cognitive 0.75
## 2 40-10_punishment overt 1
## 3 40-10_reward cognitive 1
## 4 40-10_reward overt 1
## 5 90-10_punishment cognitive 1
## 6 90-10_punishment overt 1
## 7 90-10_reward cognitive 0
## 8 90-10_reward overt 1Summarize the correct and incorrect q-values
corr_q_vals_summ <-
s1_train_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_correct_Q_vals))## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.corr_q_vals_summ$valence <- factor(corr_q_vals_summ$valence, levels=c("reward", "punishment"))
corr_q_vals_summ[corr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"incorr_q_vals_summ <-
s1_train_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_incorrect_Q_vals))## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.incorr_q_vals_summ$valence <- factor(incorr_q_vals_summ$valence, levels=c("reward", "punishment"))
incorr_q_vals_summ[incorr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"Summarize the simulated performance during learning
sim_perf_summ <-
s1_train_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_corrs))## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.sim_perf_summ$valence <- factor(sim_perf_summ$valence, levels=c("reward", "punishment"))
sim_perf_summ[sim_perf_summ$cond=="cog", "cond"] <- "cognitive"Take the difference in q-values
assert(all(corr_q_vals_summ$trial_within_condition==incorr_q_vals_summ$trial_within_condition))
assert(all(corr_q_vals_summ$valence==incorr_q_vals_summ$valence))
assert(all(corr_q_vals_summ$probability==incorr_q_vals_summ$probability))
corr_q_vals_summ$diff <- corr_q_vals_summ$m-incorr_q_vals_summ$ma <-
ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=m, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.9, color="blue", size=1.5, linetype="longdash") +
geom_hline(yintercept=-.1, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Q-value") +
ggtitle("Correct") +
xlab("Stim iteration") + ylim(-1, 1)b <-
ggplot(incorr_q_vals_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=m, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.1, color="blue", size=1.5, linetype="longdash") +
geom_hline(yintercept=-.9, color="red", size=1.5, linetype="longdash") +
ga + ap + tp + ft + tol +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ggtitle("Incorrect") +
ylab("") +
xlab("Stim iteration") + ylim(-1, 1)c <-
ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=diff, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Difference") +
xlab("Stim iteration") + ylim(-1, 1)
# For legend
# ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),
# aes(x=trial_within_condition, y=diff, color=valence)) +
# geom_hline(yintercept=0, color="black", size=2) +
# geom_line(size=3) + facet_wrap( ~ cond) +
# scale_color_manual(values=c("blue", "red")) +
# geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
# ga + ap + ft +
# theme(legend.text = element_text(size = 30),
# legend.title = element_blank(),
# legend.key.size = unit(2.5, 'lines')) +
# theme(plot.title = element_text(size = 30, hjust = .5)) +
# ylab("") +
# ggtitle("Difference") +
# xlab("Stim iteration") + ylim(-1, 1)d <- ggplot(sim_perf_summ %>% filter(probability=="90-10"),
aes(x=trial_within_condition, y=m, color=valence)) +
geom_hline(yintercept=0.5, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ cond) +
scale_color_manual(values=c("blue", "red")) +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Proportion correct") +
ggtitle("Learning") +
xlab("Stim iteration") sim_test_perf_summ <-
s1_test_sim_m1_eb %>% filter(ID==25) %>%
group_by(cond, probability, valence) %>% summarize(m=mean(sit_sim_corrs))## `summarise()` has grouped output by 'cond', 'probability'. You can override
## using the `.groups` argument.sim_test_perf_summ$valence <- factor(sim_test_perf_summ$valence, levels=c("reward", "punishment"))
sim_test_perf_summ[sim_test_perf_summ$cond=="cog", "cond"] <- "cognitive"#sim_test_perf_summ %>% filter(probability=="90-10")
e <- ggplot(sim_test_perf_summ %>% filter(probability=="90-10"),
aes(x=valence, y=m, color=valence)) +
geom_hline(yintercept=0.5, color="black", size=2) +
geom_bar(stat="identity", fill="white", size=3) +
facet_wrap(~cond) +
scale_color_manual(values=c("blue", "red")) +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Proportion correct") +
ggtitle("Test") +
xlab("") + theme(axis.ticks.x=element_blank(), axis.text.x=element_blank())#+ ylim(.48, 1) +q_vals_comb <- a + b + cq_vals_comb
sim_perf_comb <- d + esim_perf_comb
#ggsave("../paper/figs/fig_parts/q_vals_plot.png", q_vals_comb, height=5, width=16, dpi=700)#ggsave("../paper/figs/fig_parts/sim_perf_comb.png", sim_perf_comb, height=6, width=14, dpi=700)Supplementary Figure 7
Plotting Q-values by valence for m3 vs. m4 (m6 vs. m7) to understand how the introduction of pessimistic priors affects performance
# Pessimistic
m3_sims_s1 <-
rm(sp, "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayToPessPrior82512.csv")
# Naive
m4_sims_s1 <-
rm(sp, "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayTo0Inits36394.csv")M6: Q nothing varied in paper
pess_priors_qv_summs <- m3_sims_s1 %>% filter(probability=="90-10") %>%
group_by(valence, trial_within_condition) %>%
summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>%
mutate(diff=m_corr-m_incorr)## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.pess_priors_qv_summs$valence <- factor(pess_priors_qv_summs$valence, levels=c("reward", "punishment"))M7: Q nothing varied, decay to 0
naive_priors_qv_summs <- m4_sims_s1 %>% filter(probability=="90-10") %>%
group_by(valence, trial_within_condition) %>%
summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>%
mutate(diff=m_corr-m_incorr)## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.naive_priors_qv_summs$valence <- factor(naive_priors_qv_summs$valence, levels=c("reward", "punishment"))Break down into correct/incorrect plots
ggplot(pess_priors_qv_summs,
aes(x=trial_within_condition, y=m_corr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Correct") +
xlab("Stim iteration") + ylim(-1, 1)
ggplot(naive_priors_qv_summs,
aes(x=trial_within_condition, y=m_corr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Correct") +
xlab("Stim iteration") + ylim(-1, 1)
ggplot(pess_priors_qv_summs,
aes(x=trial_within_condition, y=m_incorr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
#geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Incorrect") +
xlab("Stim iteration") + ylim(-1, 1)
ggplot(naive_priors_qv_summs,
aes(x=trial_within_condition, y=m_incorr, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
# geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Incorrect") +
xlab("Stim iteration") + ylim(-1, 1)
pess_plot <- ggplot(pess_priors_qv_summs,
aes(x=trial_within_condition, y=diff, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("Q-value Difference: \n Correct vs. Incorrect") +
ggtitle("Pessimistic Priors") +
xlab("Stim iteration") + ylim(0, 1)naive_plot <- ggplot(naive_priors_qv_summs,
aes(x=trial_within_condition, y=diff, color=valence)) +
geom_hline(yintercept=0, color="black", size=2) +
geom_line(size=3) + facet_wrap( ~ valence) +
scale_color_manual(values=c("blue", "red")) +
geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
ga + ap + tol + ft +
theme(plot.title = element_text(size = 30, hjust = .5)) +
ylab("") +
ggtitle("Naive Priors (All Zero)") +
xlab("Stim iteration") + ylim(0, 1)simple_qv_plots_comb <- pess_plot + naive_plotsimple_qv_plots_comb 
No powerpoint — loaded directly into supplemental
#ggsave("../paper/figs/fig_parts/simple_qv_plots_comb.png", simple_qv_plots_comb, height=7, width=14, dpi=700)PST/Stimulus valuation phase analyses
Used rew history because Q-values are for specific Q(s,a)s whereas selection in the PST phase is between stimuli ie. Q(s, :)
# Take just the trials within condition
within_cond_pst_s1 <- s1_pst %>% filter(test_condition == "cogn_cogn" | test_condition == "overt_overt")
within_cond_pst_s1$test_condition <- factor(within_cond_pst_s1$test_condition, levels=c("overt_overt", "cogn_cogn"))
contrasts(within_cond_pst_s1$test_condition) <- c(-.5, .5)
head(within_cond_pst_s1$test_condition)## [1] overt_overt overt_overt overt_overt overt_overt overt_overt overt_overt
## attr(,"contrasts")
## [,1]
## overt_overt -0.5
## cogn_cogn 0.5
## Levels: overt_overt cogn_cogn# Both singular
# summary(m0_s1 <- glmer(resp_num ~ scale(left_min_right) + ( scale(left_min_right)|ID),
# data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
#
# summary(m0_s1 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (scale(left_min_right) |ID),
# data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))Using model without random intercepts because not singular
summary(m0_s1 <- glmer(resp_num ~ scale(left_min_right) + (0 + scale(left_min_right)|ID),
data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 6456.5 6476.6 -3225.2 6450.5 5983
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3395 -0.7468 0.1149 0.7468 4.5130
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.538 1.24
## Number of obs: 5986, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11210 0.03073 3.648 0.000264 ***
## scale(left_min_right) 1.41202 0.12083 11.686 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scl(lft_m_) 0.015summary(m1_s1 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 6431.2 6464.7 -3210.6 6421.2 5981
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.1380 -0.7457 0.1136 0.7448 4.1139
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.555 1.247
## Number of obs: 5986, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11079 0.03085 3.592 0.000328 ***
## scale(left_min_right) 1.42529 0.12149 11.731 < 2e-16 ***
## test_condition1 0.09218 0.06167 1.495 0.134952
## scale(left_min_right):test_condition1 -0.36515 0.07029 -5.195 2.05e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sc(__) tst_c1
## scl(lft_m_) 0.015
## test_cndtn1 -0.041 0.005
## scl(l__):_1 0.003 -0.035 0.025# Take just the trials within condition
within_cond_pst_s2 <- s2_pst %>% filter(test_condition == "cognitive_cognitive" | test_condition == "overt_overt")
within_cond_pst_s2$test_condition <-
factor(within_cond_pst_s2$test_condition, levels=c("overt_overt", "cognitive_cognitive"))
contrasts(within_cond_pst_s2$test_condition) <- c(-.5, .5)
head(within_cond_pst_s2$test_condition)## [1] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## [4] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## attr(,"contrasts")
## [,1]
## overt_overt -0.5
## cognitive_cognitive 0.5
## Levels: overt_overt cognitive_cognitiveUse the comparable model
summary(m0_s2 <- glmer(resp_num ~ scale(left_min_right) + (0 + scale(left_min_right)|ID),
data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7161.3 7181.7 -3577.7 7155.3 6616
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.1402 -0.7899 0.1255 0.7803 4.7614
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.507 1.228
## Number of obs: 6619, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11126 0.02914 3.818 0.000134 ***
## scale(left_min_right) 1.35153 0.11346 11.912 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scl(lft_m_) 0.014summary(m1_s2 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |
## ID)
## Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7163.5 7197.5 -3576.8 7153.5 6614
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.0975 -0.7881 0.1258 0.7749 4.9593
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID scale(left_min_right) 1.509 1.228
## Number of obs: 6619, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11101 0.02914 3.809 0.00014 ***
## scale(left_min_right) 1.35180 0.11350 11.910 < 2e-16 ***
## test_condition1 0.06205 0.05826 1.065 0.28689
## scale(left_min_right):test_condition1 -0.05321 0.06672 -0.797 0.42519
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sc(__) tst_c1
## scl(lft_m_) 0.014
## test_cndtn1 -0.008 0.004
## scl(l__):_1 0.006 -0.002 0.027# Singular
# summary(m1_s2 <- glmer(resp_num ~ scale(left_min_right)*test_condition + (0 + scale(left_min_right)*test_condition |ID),
# data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))p_s1 <-
sjPlot::plot_model(m1_s1, type="pred", terms = c("left_min_right [all]", "test_condition"))
p_fin_1 <- p_s1 + ap + tp + ga +
xlab("Left minus right reward history") + ylab("") + ggtitle("") +
theme(plot.title = element_text(hjust=0)) +
ylab("Chance of picking left") +
xlab("Left minus right reward history") +
ggtitle("Predicted probabilities") +
scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) + tol## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.p_s2 <-
sjPlot::plot_model(m1_s2, type="pred", terms = c("left_min_right [all]", "test_condition"))
p_fin_s2 <- p_s2 + ap + tp + ga +
xlab("Left minus right reward history") + ylab("") + ggtitle("") +
theme(plot.title = element_text(hjust=0)) +
ylab("Chance of picking left") +
xlab("Left minus right reward history") +
#ggtitle("Predicted probabilities") +
scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) +
theme(legend.text = element_text(size = 18),
legend.title = element_blank(),
legend.key.size = unit(1.3, 'lines'))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.p_fin_1 + p_fin_s2
ranef(m1_s2)## $ID
## scale(left_min_right)
## 1 -7.638427e-02
## 2 -1.911542e-01
## 3 -1.535659e+00
## 4 1.398870e-01
## 5 1.763300e+00
## 6 -1.465558e+00
## 7 -5.962328e-01
## 8 8.459447e-01
## 9 -2.737737e-01
## 10 8.023781e-01
## 11 -1.789991e+00
## 12 -3.571217e-01
## 13 -1.554381e+00
## 14 6.793301e-01
## 15 3.156066e-01
## 16 -1.700867e-01
## 17 -1.283673e+00
## 18 -9.527213e-01
## 19 -7.157713e-01
## 20 -4.166149e-01
## 21 -1.332563e+00
## 22 -9.936151e-01
## 23 4.158967e-01
## 24 -1.262534e+00
## 25 -8.249032e-01
## 26 1.347361e+00
## 27 -1.197001e+00
## 28 1.846921e+00
## 29 -1.252705e+00
## 30 6.828607e-01
## 31 -1.364753e+00
## 32 -1.792755e+00
## 33 8.582025e-01
## 34 -1.775080e-01
## 35 -5.372453e-02
## 36 6.800819e-01
## 37 1.997015e-01
## 38 -5.007704e-01
## 39 1.040108e+00
## 40 2.002604e+00
## 41 -1.831436e+00
## 42 6.286168e-01
## 43 1.583285e+00
## 44 -7.409308e-01
## 45 1.068599e+00
## 46 -3.780431e-01
## 47 6.087241e-01
## 48 2.624062e-01
## 49 1.076875e+00
## 50 -1.681045e+00
## 51 7.740095e-02
## 52 1.276172e-01
## 53 -8.255441e-01
## 54 1.501652e+00
## 55 -1.198909e+00
## 56 1.231464e+00
## 57 -7.374282e-01
## 58 -1.044010e+00
## 59 1.330937e+00
## 60 9.092997e-01
## 61 6.311263e-03
## 62 1.141079e+00
## 63 -6.813984e-01
## 64 -6.699732e-01
## 65 -1.327898e+00
## 66 -1.787861e-01
## 67 -1.396548e-01
## 68 -3.535874e-01
## 69 1.867799e-01
## 70 -1.221560e+00
## 71 -5.083901e-01
## 72 9.698690e-01
## 73 6.657423e-01
## 74 -7.371641e-01
## 75 -1.415828e+00
## 76 -1.264396e+00
## 77 1.724080e+00
## 78 -1.221658e+00
## 79 1.657571e+00
## 80 -1.628377e+00
## 81 2.921264e-01
## 82 -2.411008e-01
## 83 -9.298301e-05
## 84 1.107018e+00
## 85 -9.797200e-02
## 86 -2.256282e+00
## 87 1.880165e+00
## 88 -1.238028e+00
## 89 -2.559056e+00
## 90 -1.093334e+00
## 91 -7.953529e-01
## 92 -4.959134e-01
## 93 1.412486e+00
## 94 -9.225457e-01
## 95 1.014909e+00
## 96 -8.901824e-01
## 97 2.365885e-01
## 98 -1.210926e+00
## 99 3.047484e-02
## 100 3.144635e-01
## 101 -1.402683e+00
## 102 -8.816298e-01
## 103 1.511560e+00
## 104 -9.874832e-01
## 105 2.173403e-01
## 106 -1.052314e+00
## 107 1.724758e+00
## 108 1.839174e+00
## 109 -1.697650e+00
## 110 -9.228177e-01
## 111 1.453432e+00
## 112 8.494656e-01
## 113 2.214903e+00
## 114 -4.507024e-01
## 115 -3.275397e-01
## 116 -1.026534e+00
## 117 7.187065e-01
## 118 2.520358e-01
## 119 2.038057e+00
## 120 -1.233438e+00
## 121 1.897789e+00
## 122 -1.107078e+00
## 123 1.644766e+00
## 124 9.238578e-01
## 125 9.303096e-02
## 126 -2.417370e-01
## 127 1.273118e+00
## 128 6.328269e-01
## 129 1.196723e+00
## 130 -1.583174e-01
## 131 -1.289770e+00
## 132 -5.052138e-01
## 133 7.200495e-01
## 134 9.689824e-01
## 135 1.159983e+00
## 136 -1.270916e+00
## 137 -2.006149e-01
## 138 1.017842e+00
##
## with conditional variances for "ID"Mixed test phase
unique(s1_pst$test_condition)## [1] "overt_cogn" "overt_overt" "cogn_cogn"mix_s1 <- s1_pst %>% filter(test_condition == "overt_cogn")# If cognitive is on the left and they chose left then they chose cognitive
mix_s1[mix_s1$left_training_cond=="cognTraining", "chose_cognitive"] <-
if_else(mix_s1[mix_s1$left_training_cond=="cognTraining", "response"]=="left", 1, 0)
# If cognitive is on the right and they chose right then they chose cognitive
mix_s1[mix_s1$right_training_cond=="cognTraining", "chose_cognitive"] <-
if_else(mix_s1[mix_s1$right_training_cond=="cognTraining", "response"]=="right", 1, 0)table(mix_s1$chose_cognitive)[2]/sum(table(mix_s1$chose_cognitive))## 1
## 0.4751506unique(s2_pst$test_condition)## [1] "overt_cognitive" "cognitive_cognitive" "overt_overt"mix_s2 <- s2_pst %>% filter(test_condition == "overt_cognitive")# If cognitive is on the left and they chose left then they chose cognitive
mix_s2[mix_s2$left_training_cond=="cognitive", "chose_cognitive"] <-
if_else(mix_s2[mix_s2$left_training_cond=="cognitive", "response"]=="left", 1, 0)
# If cognitive is on the right and they chose right then they chose cognitive
mix_s2[mix_s2$right_training_cond=="cognitive", "chose_cognitive"] <-
if_else(mix_s2[mix_s2$right_training_cond=="cognitive", "response"]=="right", 1, 0)table(mix_s2$chose_cognitive)[2]/sum(table(mix_s2$chose_cognitive))## 1
## 0.4231989mixed_summs_s1 <- mix_s1 %>% group_by(choice_type_plotting) %>%
summarize(m=mean(chose_cognitive), n())
mixed_summs_s1## # A tibble: 4 × 3
## choice_type_plotting m `n()`
## <chr> <dbl> <int>
## 1 "freq-P\nfreq-P" 0.490 997
## 2 "freq-R\nfreq-R" 0.428 997
## 3 "infreq-P\ninfreq-P" 0.553 992
## 4 "infreq-R\ninfreq-R" 0.429 998mixed_summs_s2 <- mix_s2 %>% group_by(choice_type_plotting) %>%
summarize(m=mean(chose_cognitive), n())
mixed_summs_s2## # A tibble: 4 × 3
## choice_type_plotting m `n()`
## <chr> <dbl> <int>
## 1 "freq-P\nfreq-P" 0.408 1104
## 2 "freq-R\nfreq-R" 0.393 1104
## 3 "infreq-P\ninfreq-P" 0.447 1102
## 4 "infreq-R\ninfreq-R" 0.445 1104summary(choice_effect_s1 <- glmer(chose_cognitive ~ 1 + ( 1 |ID),
data=mix_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
## Data: mix_s1
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 5269.0 5281.6 -2632.5 5265.0 3982
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4732 -0.8937 -0.4973 0.9600 2.3488
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.5492 0.7411
## Number of obs: 3984, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.10926 0.07443 -1.468 0.142summary(choice_effect_s2 <- glmer(chose_cognitive ~ 1 + ( 1|ID),
data=mix_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
## Data: mix_s2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 5829.4 5842.2 -2912.7 5825.4 4412
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5760 -0.7917 -0.5754 1.0417 2.2498
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.3916 0.6258
## Number of obs: 4414, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.33769 0.06215 -5.433 5.53e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Questionnaires
s1_qs <- read.csv("../data/cleaned_files/s1_qs_deidentified.csv")
s2_qs <- read.csv("../data/cleaned_files/s2_qs_deidentified.csv")s1_demogs <- read.csv("../data/cleaned_files/s1_demogs_deIDed.csv")
s2_demogs <- read.csv("../data/cleaned_files/s2_demogs_deIDed.csv")# names(s1_qs)
# names(s2_qs)Exclude pts who answered both catch questions incorrectly
S1 all correct so no exclusions
table(s1_qs$Q6_20) # All correct##
## A little bit
## 123table(s1_qs$Q7_19) # All correct##
## Very false for me
## 1234/4 — now removing the pt who skipped all items here (leaving n=122 completers). This mayslightly influenced PTQ and RRS imputations; statistic now corrected to p’s > .07 (originally .10)
s1_qs <- s1_qs[-66, ]table(s2_qs$Q6_20) # All correct##
## A little bit
## 137table(s2_qs$Q7_19) # One pt incorrect but they answered the other catch item correctly, so retained (given recs in lit that excluding on a single catch item can bias results)##
## Very false for me Very true for me
## 136 1#s2_qs %>% filter(Q7_19 == "Very true for me")From notes:
Q2 is the PTQ
Q6 is the MASQ — 23 items, catch is 20 (confirmed on Qualtrics) - but
permission issue explained in revision hence not reported
Q7 is the BIS/BAS — 25 items, catch is 19 (confirmed on Qualtrics)
- some are not part of the scale
Q8 is the RRS — 10 items, no catch
S1 questionnaires
Get and recode PTQ and RRS
s1_ptq <- s1_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
#s1_ptq$ID <- s1_qs$ID 4/4 removing to make sure not in imputation
s1_ptq_copy <- s1_ptq
s1_ptq_copy$ID <- s1_qs$ID
s1_ptq[s1_ptq=="Never"] <- 0
s1_ptq[s1_ptq=="Rarely"] <- 1
s1_ptq[s1_ptq=="Sometimes"] <- 2
s1_ptq[s1_ptq=="Often"] <- 3
s1_ptq[s1_ptq=="Almost always"] <- 4s1_rrs <-
s1_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
#s1_rrs$ID <- s1_qs$ID 4/4 removing to make sure not in imputation
s1_rrs_copy <- s1_rrs
s1_rrs_copy$ID <- s1_qs$ID
s1_rrs[s1_rrs=="Almost never"] <- 1
s1_rrs[s1_rrs=="Sometimes"] <- 2
s1_rrs[s1_rrs=="Often"] <- 3
s1_rrs[s1_rrs=="Almost always"] <- 4For revision: BIS/BAS
Exclude catch but keep filler items for now
#s1_qs %>% select(contains("Q7_")) %>% select("Q7_19") # confirmed this is catch
s1_bis_bas <- s1_qs %>% select(contains("Q7_")) %>% select(-"Q7_19") “Items other than 2 and 22 are reverse coded”
The two BIS items that are forward coded (1= very true, 4=very false) — Q7_23 corresponds to item 22 in actual scale because of the catch item
fwd_code_bis_bas <- s1_bis_bas %>% select("Q7_2", "Q7_23")
fwd_code_bis_bas[fwd_code_bis_bas=="Very true for me"] <- 1
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat true for me"] <- 2
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat false for me"] <- 3
fwd_code_bis_bas[fwd_code_bis_bas=="Very false for me"] <- 4Sanity check that positively correlate
cor.test(as.numeric(fwd_code_bis_bas$Q7_2), as.numeric(fwd_code_bis_bas$Q7_23))##
## Pearson's product-moment correlation
##
## data: as.numeric(fwd_code_bis_bas$Q7_2) and as.numeric(fwd_code_bis_bas$Q7_23)
## t = 8.8719, df = 120, p-value = 8.218e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5085008 0.7259186
## sample estimates:
## cor
## 0.6293718The rest are reverse coded
rev_code_bis_bas <- s1_bis_bas %>% select(-c("Q7_2", "Q7_23"))
rev_code_bis_bas[rev_code_bis_bas=="Very true for me"] <- 4
rev_code_bis_bas[rev_code_bis_bas=="Somewhat true for me"] <- 3
rev_code_bis_bas[rev_code_bis_bas=="Somewhat false for me"] <- 2
rev_code_bis_bas[rev_code_bis_bas=="Very false for me"] <- 1s1_bis_bas <- data.frame(fwd_code_bis_bas, rev_code_bis_bas)s1_bis_bas_num <- foreach (i = 1:ncol(s1_bis_bas)) %do% {
data.frame(as.numeric(unlist(s1_bis_bas[i]))) %>% setNames(names(s1_bis_bas[i]))
}%>% bind_cols()Sanity check BIS fwd and reverse code items postiively correlate
cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_23)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_23
## t = 8.8719, df = 120, p-value = 8.218e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5085008 0.7259186
## sample estimates:
## cor
## 0.6293718cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_8)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_8
## t = 6.3919, df = 120, p-value = 3.264e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3583137 0.6256823
## sample estimates:
## cor
## 0.5039744cor.test(s1_bis_bas_num$Q7_23, s1_bis_bas_num$Q7_8)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_23 and s1_bis_bas_num$Q7_8
## t = 4.8078, df = 120, p-value = 4.47e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2413693 0.5409989
## sample estimates:
## cor
## 0.4018868# Item 24 because shifted by 1 due to catch
cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_25)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_25
## t = 8.6941, df = 120, p-value = 2.152e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4990482 0.7198710
## sample estimates:
## cor
## 0.6216608cor.test(s1_bis_bas_num$Q7_23, s1_bis_bas_num$Q7_25)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_23 and s1_bis_bas_num$Q7_25
## t = 5.8332, df = 120, p-value = 4.707e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3188964 0.5978260
## sample estimates:
## cor
## 0.4700135cor.test(s1_bis_bas_num$Q7_2, s1_bis_bas_num$Q7_16)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_2 and s1_bis_bas_num$Q7_16
## t = 6.2646, df = 120, p-value = 6.056e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3495145 0.6195225
## sample estimates:
## cor
## 0.4964321cor.test(s1_bis_bas_num$Q7_23, s1_bis_bas_num$Q7_16)##
## Pearson's product-moment correlation
##
## data: s1_bis_bas_num$Q7_23 and s1_bis_bas_num$Q7_16
## t = 6.0606, df = 120, p-value = 1.611e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3351901 0.6094231
## sample estimates:
## cor
## 0.4841061Now remove filler items — all before catch so all items correspond to actual num in scale
s1_bis_bas <- s1_bis_bas_num %>% select(-c(
"Q7_1",
"Q7_6",
"Q7_11",
"Q7_17")
)BIS
s1_bis <- s1_bis_bas %>% select(
# Before catch item
c("Q7_2", "Q7_8", "Q7_13", "Q7_16"),
# After catch so one shifted
c("Q7_20", "Q7_23", "Q7_25"))Just 1 missing so impute
length(which(is.na(s1_bis[, 1:7])))## [1] 1length(which(is.na(s1_bis[, 1:7])))/(dim(s1_bis[, 1:7])[1]*dim(s1_bis[, 1:7])[2]) ## [1] 0.00117096s1_bis_final <- foreach (i = 1:nrow(s1_bis)) %do% {
if (any(is.na(s1_bis[i, ]))) {
s1_bis[i, is.na(s1_bis[i, ])] <- mean(unlist(s1_bis[i, ]), na.rm=TRUE)
}
s1_bis[i, ]
}%>% bind_rows()Sanity check all positively correlated
cor(s1_bis_final)## Q7_2 Q7_8 Q7_13 Q7_16 Q7_20 Q7_23 Q7_25
## Q7_2 1.0000000 0.5039744 0.5156167 0.4964321 0.5391917 0.6293718 0.6216608
## Q7_8 0.5039744 1.0000000 0.6775552 0.5204724 0.5118818 0.4018868 0.5865127
## Q7_13 0.5156167 0.6775552 1.0000000 0.5338191 0.5322763 0.4445993 0.4822497
## Q7_16 0.4964321 0.5204724 0.5338191 1.0000000 0.5482658 0.4841061 0.5369988
## Q7_20 0.5391917 0.5118818 0.5322763 0.5482658 1.0000000 0.3565339 0.6321968
## Q7_23 0.6293718 0.4018868 0.4445993 0.4841061 0.3565339 1.0000000 0.4700135
## Q7_25 0.6216608 0.5865127 0.4822497 0.5369988 0.6321968 0.4700135 1.0000000And high alpha
psych::alpha(s1_bis_final)##
## Reliability analysis
## Call: psych::alpha(x = s1_bis_final)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.88 0.89 0.88 0.53 7.7 0.016 3 0.7 0.52
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.85 0.88 0.91
## Duhachek 0.85 0.88 0.91
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q7_2 0.86 0.86 0.86 0.51 6.4 0.020 0.0067 0.52
## Q7_8 0.86 0.87 0.86 0.52 6.5 0.019 0.0053 0.53
## Q7_13 0.86 0.87 0.86 0.52 6.6 0.019 0.0061 0.52
## Q7_16 0.87 0.87 0.87 0.53 6.7 0.019 0.0082 0.52
## Q7_20 0.87 0.87 0.86 0.53 6.7 0.019 0.0055 0.52
## Q7_23 0.88 0.88 0.87 0.55 7.3 0.017 0.0031 0.53
## Q7_25 0.86 0.86 0.86 0.51 6.3 0.019 0.0062 0.52
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q7_2 122 0.80 0.80 0.76 0.71 3.0 0.90
## Q7_8 122 0.79 0.78 0.74 0.69 2.8 1.01
## Q7_13 122 0.79 0.78 0.74 0.69 2.9 0.99
## Q7_16 122 0.76 0.76 0.71 0.67 2.8 0.88
## Q7_20 122 0.75 0.76 0.72 0.66 3.3 0.81
## Q7_23 122 0.71 0.70 0.64 0.59 2.8 1.00
## Q7_25 122 0.79 0.80 0.77 0.72 3.2 0.83
##
## Non missing response frequency for each item
## 1 2 2.66666666666667 3 4 miss
## Q7_2 0.07 0.18 0.00 0.39 0.36 0
## Q7_8 0.13 0.25 0.00 0.33 0.29 0
## Q7_13 0.11 0.20 0.00 0.34 0.34 0
## Q7_16 0.07 0.30 0.00 0.41 0.22 0
## Q7_20 0.06 0.06 0.01 0.45 0.43 0
## Q7_23 0.12 0.25 0.00 0.34 0.29 0
## Q7_25 0.04 0.14 0.00 0.41 0.41 0s1_bis_final$bis_sum <- rowSums(s1_bis_final)
s1_bis_final$ID <- s1_qs$IDThe remaining items are BAS
s1_bas <- s1_bis_bas %>% select(
-c("Q7_2", "Q7_8", "Q7_13", "Q7_16",
"Q7_20", "Q7_23", "Q7_25"))Just 3 missing so impute
length(which(is.na(s1_bas[, 1:ncol(s1_bas)])))## [1] 3length(which(is.na(s1_bas[, 1:ncol(s1_bas)])))/(dim(s1_bas[, 1:ncol(s1_bas)])[1]*dim(s1_bas[, 1:ncol(s1_bas)])[2]) ## [1] 0.001891551s1_bas_final <- foreach (i = 1:nrow(s1_bas)) %do% {
if (any(is.na(s1_bas[i, ]))) {
s1_bas[i, is.na(s1_bas[i, ])] <- mean(unlist(s1_bas[i, ]), na.rm=TRUE)
}
s1_bas[i, ]
}%>% bind_rows()Sanity check all positively correlated
cor(s1_bas_final)## Q7_3 Q7_4 Q7_5 Q7_7 Q7_9 Q7_10 Q7_12
## Q7_3 1.0000000 0.3505651 0.5194718 0.3841874 0.6711250 0.5025063 0.6169502
## Q7_4 0.3505651 1.0000000 0.2879864 0.3390419 0.2583916 0.1969002 0.3028181
## Q7_5 0.5194718 0.2879864 1.0000000 0.4144983 0.4738954 0.6536637 0.5278935
## Q7_7 0.3841874 0.3390419 0.4144983 1.0000000 0.3083571 0.3054471 0.2941940
## Q7_9 0.6711250 0.2583916 0.4738954 0.3083571 1.0000000 0.4900888 0.7280516
## Q7_10 0.5025063 0.1969002 0.6536637 0.3054471 0.4900888 1.0000000 0.3838347
## Q7_12 0.6169502 0.3028181 0.5278935 0.2941940 0.7280516 0.3838347 1.0000000
## Q7_14 0.5878563 0.3752357 0.5274854 0.6148655 0.5484696 0.3908887 0.6012681
## Q7_15 0.2808272 0.1617566 0.4120492 0.2799688 0.3525170 0.5021741 0.3237509
## Q7_18 0.2918551 0.1745012 0.1677867 0.5125142 0.2076797 0.1884286 0.2343946
## Q7_21 0.5414331 0.1512069 0.5872976 0.3294216 0.4492972 0.5246798 0.3926814
## Q7_22 0.5064174 0.2087005 0.3996587 0.1100529 0.6215593 0.3941751 0.5789322
## Q7_24 0.3864836 0.2818568 0.3438651 0.5351943 0.2465222 0.1690507 0.3441865
## Q7_14 Q7_15 Q7_18 Q7_21 Q7_22 Q7_24
## Q7_3 0.5878563 0.2808272 0.2918551 0.5414331 0.5064174 0.3864836
## Q7_4 0.3752357 0.1617566 0.1745012 0.1512069 0.2087005 0.2818568
## Q7_5 0.5274854 0.4120492 0.1677867 0.5872976 0.3996587 0.3438651
## Q7_7 0.6148655 0.2799688 0.5125142 0.3294216 0.1100529 0.5351943
## Q7_9 0.5484696 0.3525170 0.2076797 0.4492972 0.6215593 0.2465222
## Q7_10 0.3908887 0.5021741 0.1884286 0.5246798 0.3941751 0.1690507
## Q7_12 0.6012681 0.3237509 0.2343946 0.3926814 0.5789322 0.3441865
## Q7_14 1.0000000 0.3246495 0.4751540 0.4521666 0.4152103 0.5021592
## Q7_15 0.3246495 1.0000000 0.2702716 0.5064189 0.3647177 0.1475499
## Q7_18 0.4751540 0.2702716 1.0000000 0.1857702 0.2346360 0.4633299
## Q7_21 0.4521666 0.5064189 0.1857702 1.0000000 0.4651129 0.3173153
## Q7_22 0.4152103 0.3647177 0.2346360 0.4651129 1.0000000 0.1378870
## Q7_24 0.5021592 0.1475499 0.4633299 0.3173153 0.1378870 1.0000000And high alpha
psych::alpha(s1_bas_final) ##
## Reliability analysis
## Call: psych::alpha(x = s1_bas_final)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.89 0.89 0.92 0.39 8.2 0.014 2.9 0.52 0.38
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.86 0.89 0.92
## Duhachek 0.86 0.89 0.92
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q7_3 0.88 0.88 0.90 0.37 7.1 0.016 0.022 0.36
## Q7_4 0.89 0.89 0.92 0.41 8.3 0.015 0.021 0.41
## Q7_5 0.88 0.88 0.90 0.38 7.2 0.016 0.022 0.36
## Q7_7 0.89 0.88 0.91 0.39 7.6 0.015 0.023 0.39
## Q7_9 0.88 0.88 0.90 0.38 7.2 0.016 0.020 0.38
## Q7_10 0.88 0.88 0.91 0.38 7.5 0.015 0.022 0.37
## Q7_12 0.88 0.88 0.90 0.38 7.2 0.016 0.021 0.38
## Q7_14 0.88 0.87 0.90 0.37 7.0 0.016 0.022 0.35
## Q7_15 0.89 0.89 0.91 0.40 7.9 0.015 0.023 0.39
## Q7_18 0.89 0.89 0.91 0.40 8.2 0.014 0.021 0.39
## Q7_21 0.88 0.88 0.91 0.38 7.4 0.016 0.023 0.37
## Q7_22 0.88 0.88 0.91 0.39 7.6 0.015 0.021 0.38
## Q7_24 0.89 0.89 0.91 0.40 7.9 0.015 0.022 0.39
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q7_3 122 0.78 0.78 0.76 0.72 2.9 0.81
## Q7_4 122 0.44 0.48 0.40 0.37 3.5 0.53
## Q7_5 122 0.75 0.74 0.72 0.68 3.0 0.83
## Q7_7 122 0.62 0.63 0.61 0.54 3.3 0.74
## Q7_9 122 0.75 0.74 0.74 0.69 2.8 0.79
## Q7_10 122 0.69 0.67 0.64 0.60 2.7 0.89
## Q7_12 122 0.74 0.74 0.73 0.68 2.8 0.76
## Q7_14 122 0.78 0.80 0.79 0.74 3.1 0.71
## Q7_15 122 0.59 0.58 0.53 0.50 2.2 0.82
## Q7_18 122 0.51 0.51 0.47 0.41 3.1 0.81
## Q7_21 122 0.71 0.69 0.67 0.63 2.6 0.85
## Q7_22 122 0.65 0.64 0.61 0.57 2.3 0.86
## Q7_24 122 0.55 0.57 0.53 0.47 3.4 0.70
##
## Non missing response frequency for each item
## 1 2 2.91666666666667 3 3.41666666666667 3.58333333333333 4
## Q7_3 0.07 0.20 0.00 0.52 0.00 0.00 0.20
## Q7_4 0.00 0.02 0.00 0.47 0.01 0.00 0.51
## Q7_5 0.04 0.20 0.00 0.43 0.00 0.00 0.32
## Q7_7 0.03 0.07 0.00 0.47 0.00 0.00 0.43
## Q7_9 0.04 0.33 0.00 0.45 0.00 0.01 0.17
## Q7_10 0.10 0.31 0.00 0.41 0.00 0.00 0.18
## Q7_12 0.03 0.29 0.00 0.49 0.00 0.00 0.19
## Q7_14 0.02 0.15 0.00 0.54 0.00 0.00 0.30
## Q7_15 0.19 0.45 0.00 0.30 0.00 0.00 0.06
## Q7_18 0.05 0.14 0.01 0.48 0.00 0.00 0.32
## Q7_21 0.09 0.33 0.00 0.43 0.00 0.00 0.16
## Q7_22 0.18 0.41 0.00 0.33 0.00 0.00 0.08
## Q7_24 0.02 0.07 0.00 0.39 0.00 0.00 0.52
## miss
## Q7_3 0
## Q7_4 0
## Q7_5 0
## Q7_7 0
## Q7_9 0
## Q7_10 0
## Q7_12 0
## Q7_14 0
## Q7_15 0
## Q7_18 0
## Q7_21 0
## Q7_22 0
## Q7_24 0s1_bas_final$bas_sum <- rowSums(s1_bas_final)
hist(s1_bas_final$bas_sum, breaks=25)
s1_bas_final$ID <- s1_qs$IDModestly inversely related
ComparePars(s1_bas_final$bas_sum, s1_bis_final$bis_sum, use_identity_line = 0)
Prep PTQ S1
Convert to numeric
s1_ptq_num <- foreach (i = 1:ncol(s1_ptq)) %do% {
data.frame(as.numeric(unlist(s1_ptq[i]))) %>% setNames(names(s1_ptq[i]))
}%>% bind_cols()Just 4 data points (< .1% of data) missing for remaining, so mean impute these
length(which(is.na(s1_ptq_num[, 1:15])))## [1] 4length(which(is.na(s1_ptq_num[, 1:15])))/(dim(s1_ptq_num[, 1:15])[1]*dim(s1_ptq_num[, 1:15])[2]) ## [1] 0.002185792s1_ptq_final <- foreach (i = 1:nrow(s1_ptq_num)) %do% {
if (any(is.na(s1_ptq_num[i, ]))) {
s1_ptq_num[i, is.na(s1_ptq_num[i, ])] <- mean(unlist(s1_ptq_num[i, ]), na.rm=TRUE)
}
s1_ptq_num[i, ]
}%>% bind_rows()hist(rowSums(s1_ptq_final[, 1:15]), breaks=25)
s1_ptq_final$ptq_sum <- rowSums(s1_ptq_final[, 1:15])s1_ptq_final$ID <- s1_qs$IDSanity check that correlated in expected irection
ComparePars(s1_ptq_final$ptq_sum, s1_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s1_ptq_final$ptq_sum, s1_bas_final$bas_sum, use_identity_line = 0)
Spot check IDs lined up properly
# s1_ptq_final %>% filter(ID == 22)
# s1_ptq_copy %>% filter(ID == 22)
# s1_qs %>% filter(ID == 22)
# s1_ptq_final %>% filter(ID == 104)
# s1_ptq_copy %>% filter(ID == 104)
# s1_qs %>% filter(ID == 104)Prep S1 RRS
s1_rrs_num <- foreach (i = 1:ncol(s1_rrs)) %do% {
data.frame(as.numeric(unlist(s1_rrs[i]))) %>% setNames(names(s1_rrs[i]))
}%>% bind_cols()length(which(is.na(s1_rrs_num[, 1:10])))## [1] 3length(which(is.na(s1_rrs_num[, 1:10])))/(dim(s1_rrs_num[, 1:10])[1]*dim(s1_rrs_num[, 1:10])[2]) ## [1] 0.002459016s1_rrs_final <- foreach (i = 1:nrow(s1_rrs_num)) %do% {
if (any(is.na(s1_rrs_num[i, ]))) {
s1_rrs_num[i, is.na(s1_rrs_num[i, ])] <- mean(unlist(s1_rrs_num[i, ]), na.rm=TRUE)
}
s1_rrs_num[i, ]
}%>% bind_rows()hist(rowSums(s1_rrs_final[, 1:10]), breaks=25)
s1_rrs_final$rrs_sum <- rowSums(s1_rrs_final[, 1:10])s1_rrs_final$ID <- s1_qs$IDs1 Spot check
# s1_rrs_final %>% filter(ID == 106)
# s1_rrs_copy %>% filter(ID == 106)
# s1_qs %>% filter(ID == 106) # Q8 is RRSSanity check RRS and PTQ are strongly correlated
ComparePars(s1_rrs_final$rrs_sum, s1_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(s1_rrs_final$rrs_sum, s1_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s1_rrs_final$rrs_sum, s1_bas_final$bas_sum, use_identity_line = 0)
Reduce to just matching IDs
m1_s1_model_red <- m1_study1_eb %>% filter(ID %in% s1_ptq_final$ID)
m1_s1_model_red$c_min_o_phi <- m1_s1_model_red$cog_phi - m1_s1_model_red$overt_phiassert(all(m1_s1_model_red$ID==s1_ptq_final$ID))
assert(all(m1_s1_model_red$ID==s1_rrs_final$ID))PTQ
ComparePars(m1_s1_model_red$c_min_o_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s1_model_red$cog_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s1_model_red$overt_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)
Same basic pattern for RRS 4/4 - this was where p value for top statistic slightly different due to imputation error in first version involving ID, adjusted in manuscript from
ComparePars(m1_s1_model_red$c_min_o_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)
# Would expect higher decay to positively correlate w more brooding — instead small and ns neg
ComparePars(m1_s1_model_red$cog_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)
ComparePars(m1_s1_model_red$overt_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)
Do same thing with age
Complete age data so don’t need to reduce
assert(all(s1_demogs$ID==m1_study1_eb$ID))m1_study1_eb$c_min_o_phi <- m1_study1_eb$cog_phi - m1_study1_eb$overt_phiAdded for second revision — Age — surprisingly not related to particularly high decay in cog. Closest to signif is actually higher decay in overt!
ComparePars(m1_study1_eb$c_min_o_phi, s1_demogs$Age, use_identity_line = 0)
ComparePars(m1_study1_eb$cog_phi, s1_demogs$Age, use_identity_line = 0)
ComparePars(m1_study1_eb$overt_phi, s1_demogs$Age, use_identity_line = 0)
Age and task performance
s1_train_perf <- s1_train %>% group_by(ID) %>% summarize(m=mean(correct))
s1_test_perf <- s1_sit %>% group_by(ID) %>% summarize(m=mean(correct))assert(all(s1_demogs$ID==s1_train_perf$ID))
assert(all(s1_demogs$ID==s1_test_perf$ID))Related to overall performance
ComparePars(s1_train_perf$m, s1_demogs$Age, use_identity_line = 0)
ComparePars(s1_test_perf$m, s1_demogs$Age, use_identity_line = 0)
But not difference in cog vs. overt performance
s1_perf <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_perf_s1 <- as.numeric(unlist(s1_perf[s1_perf$condition=="overt","m"]-
s1_perf[s1_perf$condition=="cognitive", "m"]))
# .. and test subject-level differences in performance
s1_perf_test <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_testperf_s1 <- as.numeric(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]-
s1_perf_test[s1_perf_test$condition=="cognitive", "m"]))ComparePars(o_min_c_perf_s1, s1_demogs$Age, use_identity_line = 0)
ComparePars(o_min_c_testperf_s1, s1_demogs$Age, use_identity_line = 0)
Added for second revision — any relationships with reward vs. punishment learning and BIS/BAS?
s1_rew_train_perf <- s1_train %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_pun_train_perf <- s1_train %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_rew_train_perf$rew_min_pun <- s1_rew_train_perf$m-s1_pun_train_perf$m
s1_rew_sit_perf <- s1_sit %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_pun_sit_perf <- s1_sit %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s1_rew_sit_perf$rew_min_pun <- s1_rew_sit_perf$m-s1_pun_sit_perf$m
ComparePars(s1_rew_train_perf$rew_min_pun, s1_rew_sit_perf$rew_min_pun, use_identity_line = 0)
Reduce to just matching IDs
s1_rew_train_perf_red <- s1_rew_train_perf %>% filter(ID %in% s1_bis_final$ID)
s1_rew_sit_perf_red <- s1_rew_sit_perf %>% filter(ID %in% s1_bis_final$ID)
ComparePars(s1_rew_train_perf_red$rew_min_pun, s1_rew_sit_perf_red$rew_min_pun, use_identity_line = 0)
assert(all(s1_rew_sit_perf_red$ID==s1_bis_final$ID))
assert(all(s1_rew_sit_perf_red$ID==s1_bas_final$ID))
assert(all(s1_rew_train_perf_red$ID==s1_bis_final$ID))
assert(all(s1_rew_train_perf_red$ID==s1_bas_final$ID))Not correlated
ComparePars(s1_rew_train_perf_red$rew_min_pun, s1_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s1_rew_train_perf_red$rew_min_pun, s1_bas_final$bas_sum, use_identity_line = 0)
ComparePars(s1_rew_sit_perf_red$rew_min_pun, s1_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s1_rew_sit_perf_red$rew_min_pun, s1_bas_final$bas_sum, use_identity_line = 0)
s2 questionnaires new code
No pts with all skips here and table above shows no exclusions
Get and recode PTQ and RRS
s2_ptq <- s2_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
#s2_ptq$ID <- s2_qs$ID 4/4 removing to make sure not in imputation
s2_ptq_copy <- s2_ptq
s2_ptq_copy$ID <- s2_qs$ID
s2_ptq[s2_ptq=="Never"] <- 0
s2_ptq[s2_ptq=="Rarely"] <- 1
s2_ptq[s2_ptq=="Sometimes"] <- 2
s2_ptq[s2_ptq=="Often"] <- 3
s2_ptq[s2_ptq=="Almost always"] <- 4s2_rrs <-
s2_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
#s2_rrs$ID <- s2_qs$ID 4/4 removing to make sure not in imputation
s2_rrs_copy <- s2_rrs
s2_rrs_copy$ID <- s2_qs$ID
s2_rrs[s2_rrs=="Almost never"] <- 1
s2_rrs[s2_rrs=="Sometimes"] <- 2
s2_rrs[s2_rrs=="Often"] <- 3
s2_rrs[s2_rrs=="Almost always"] <- 4For revision: BIS/BAS
Exclude catch but keep filler items for now
#s2_qs %>% select(contains("Q7_")) %>% select("Q7_19") # confirmed this is catch
s2_bis_bas <- s2_qs %>% select(contains("Q7_")) %>% select(-"Q7_19") “Items other than 2 and 22 are reverse coded”
The two BIS items that are forward coded (1= very true, 4=very false) — Q7_23 corresponds to item 22 in actual scale because of the catch item
fwd_code_bis_bas <- s2_bis_bas %>% select("Q7_2", "Q7_23")
fwd_code_bis_bas[fwd_code_bis_bas=="Very true for me"] <- 1
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat true for me"] <- 2
fwd_code_bis_bas[fwd_code_bis_bas=="Somewhat false for me"] <- 3
fwd_code_bis_bas[fwd_code_bis_bas=="Very false for me"] <- 4Sanity check that positively correlate
cor.test(as.numeric(fwd_code_bis_bas$Q7_2), as.numeric(fwd_code_bis_bas$Q7_23))##
## Pearson's product-moment correlation
##
## data: as.numeric(fwd_code_bis_bas$Q7_2) and as.numeric(fwd_code_bis_bas$Q7_23)
## t = 8.41, df = 134, p-value = 5.365e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4654952 0.6880297
## sample estimates:
## cor
## 0.5877714The rest are reverse coded
rev_code_bis_bas <- s2_bis_bas %>% select(-c("Q7_2", "Q7_23"))
rev_code_bis_bas[rev_code_bis_bas=="Very true for me"] <- 4
rev_code_bis_bas[rev_code_bis_bas=="Somewhat true for me"] <- 3
rev_code_bis_bas[rev_code_bis_bas=="Somewhat false for me"] <- 2
rev_code_bis_bas[rev_code_bis_bas=="Very false for me"] <- 1s2_bis_bas <- data.frame(fwd_code_bis_bas, rev_code_bis_bas)s2_bis_bas_num <- foreach (i = 1:ncol(s2_bis_bas)) %do% {
data.frame(as.numeric(unlist(s2_bis_bas[i]))) %>% setNames(names(s2_bis_bas[i]))
}%>% bind_cols()Sanity check BIS fwd and reverse code items positively correlate
cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_23)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_23
## t = 8.41, df = 134, p-value = 5.365e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4654952 0.6880297
## sample estimates:
## cor
## 0.5877714cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_8)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_8
## t = 9.1229, df = 134, p-value = 9.631e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5030616 0.7130186
## sample estimates:
## cor
## 0.6189783cor.test(s2_bis_bas_num$Q7_23, s2_bis_bas_num$Q7_8)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_23 and s2_bis_bas_num$Q7_8
## t = 6.2383, df = 135, p-value = 5.312e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3316317 0.5936545
## sample estimates:
## cor
## 0.473037# Item 24 because shifted by 1 due to catch
cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_25)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_25
## t = 10.117, df = 134, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5507643 0.7439995
## sample estimates:
## cor
## 0.6580849cor.test(s2_bis_bas_num$Q7_23, s2_bis_bas_num$Q7_25)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_23 and s2_bis_bas_num$Q7_25
## t = 6.5828, df = 135, p-value = 9.4e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3545357 0.6102077
## sample estimates:
## cor
## 0.4929404cor.test(s2_bis_bas_num$Q7_2, s2_bis_bas_num$Q7_16)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_2 and s2_bis_bas_num$Q7_16
## t = 12.059, df = 134, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6295342 0.7934011
## sample estimates:
## cor
## 0.7214175cor.test(s2_bis_bas_num$Q7_23, s2_bis_bas_num$Q7_16)##
## Pearson's product-moment correlation
##
## data: s2_bis_bas_num$Q7_23 and s2_bis_bas_num$Q7_16
## t = 8.874, df = 135, p-value = 3.773e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4890406 0.7031120
## sample estimates:
## cor
## 0.6069724Now remove filler items — all before catch so all items corresopnd to actual num in scale
s2_bis_bas <- s2_bis_bas_num %>% select(-c(
"Q7_1",
"Q7_6",
"Q7_11",
"Q7_17")
)BIS
s2_bis <- s2_bis_bas %>% select(
# Before catch item
c("Q7_2", "Q7_8", "Q7_13", "Q7_16"),
# After catch so one shifted
c("Q7_20", "Q7_23", "Q7_25"))Just 1 missing so impute
length(which(is.na(s2_bis[, 1:7])))## [1] 1length(which(is.na(s2_bis[, 1:7])))/(dim(s2_bis[, 1:7])[1]*dim(s2_bis[, 1:7])[2]) ## [1] 0.001042753s2_bis_final <- foreach (i = 1:nrow(s2_bis)) %do% {
if (any(is.na(s2_bis[i, ]))) {
s2_bis[i, is.na(s2_bis[i, ])] <- mean(unlist(s2_bis[i, ]), na.rm=TRUE)
}
s2_bis[i, ]
}%>% bind_rows()Sanity check all positively correlated
cor(s2_bis_final)## Q7_2 Q7_8 Q7_13 Q7_16 Q7_20 Q7_23 Q7_25
## Q7_2 1.0000000 0.6189294 0.5450794 0.7215200 0.5609502 0.5875614 0.6590967
## Q7_8 0.6189294 1.0000000 0.6954469 0.6496518 0.6489944 0.4730370 0.7506610
## Q7_13 0.5450794 0.6954469 1.0000000 0.6336202 0.7315979 0.3881398 0.6686291
## Q7_16 0.7215200 0.6496518 0.6336202 1.0000000 0.6550656 0.6069724 0.7005364
## Q7_20 0.5609502 0.6489944 0.7315979 0.6550656 1.0000000 0.3751194 0.7093946
## Q7_23 0.5875614 0.4730370 0.3881398 0.6069724 0.3751194 1.0000000 0.4929404
## Q7_25 0.6590967 0.7506610 0.6686291 0.7005364 0.7093946 0.4929404 1.0000000And high alpha
psych::alpha(s2_bis_final)##
## Reliability analysis
## Call: psych::alpha(x = s2_bis_final)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.92 0.92 0.92 0.61 11 0.011 3 0.76 0.65
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.89 0.92 0.94
## Duhachek 0.90 0.92 0.94
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q7_2 0.91 0.90 0.90 0.61 9.5 0.012 0.0147 0.65
## Q7_8 0.90 0.90 0.90 0.60 9.1 0.013 0.0127 0.63
## Q7_13 0.91 0.91 0.90 0.61 9.5 0.012 0.0106 0.65
## Q7_16 0.90 0.90 0.90 0.59 8.8 0.013 0.0142 0.62
## Q7_20 0.91 0.90 0.90 0.61 9.5 0.012 0.0101 0.63
## Q7_23 0.92 0.92 0.92 0.66 11.8 0.010 0.0034 0.66
## Q7_25 0.90 0.90 0.90 0.59 8.7 0.013 0.0119 0.62
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q7_2 137 0.82 0.82 0.78 0.75 3.1 0.91
## Q7_8 137 0.85 0.85 0.82 0.78 2.8 0.99
## Q7_13 137 0.82 0.81 0.78 0.74 3.0 0.99
## Q7_16 137 0.87 0.87 0.85 0.81 2.7 0.98
## Q7_20 137 0.82 0.82 0.79 0.75 3.2 0.86
## Q7_23 137 0.68 0.69 0.61 0.58 2.9 0.85
## Q7_25 137 0.87 0.87 0.85 0.82 3.0 0.95
##
## Non missing response frequency for each item
## 1 2 3 3.66666666666667 4 miss
## Q7_2 0.06 0.18 0.33 0.01 0.43 0
## Q7_8 0.11 0.26 0.33 0.00 0.30 0
## Q7_13 0.09 0.21 0.31 0.00 0.39 0
## Q7_16 0.13 0.26 0.37 0.00 0.24 0
## Q7_20 0.05 0.12 0.38 0.00 0.45 0
## Q7_23 0.04 0.29 0.39 0.00 0.27 0
## Q7_25 0.09 0.17 0.36 0.00 0.39 0s2_bis_final$bis_sum <- rowSums(s2_bis_final)
s2_bis_final$ID <- s2_qs$IDThe remaining items are BAS
s2_bas <- s2_bis_bas %>% select(
-c("Q7_2", "Q7_8", "Q7_13", "Q7_16",
"Q7_20", "Q7_23", "Q7_25"))Nothing missing
s2_bas_final <- s2_basSanity check all positively correlated
cor(s2_bas_final)## Q7_3 Q7_4 Q7_5 Q7_7 Q7_9 Q7_10 Q7_12
## Q7_3 1.0000000 0.3185535 0.3257158 0.4426593 0.5869646 0.1445450 0.5553761
## Q7_4 0.3185535 1.0000000 0.3380587 0.5274478 0.2944064 0.1264353 0.3069515
## Q7_5 0.3257158 0.3380587 1.0000000 0.4129971 0.3148678 0.4485842 0.2683164
## Q7_7 0.4426593 0.5274478 0.4129971 1.0000000 0.3433022 0.2334689 0.3594777
## Q7_9 0.5869646 0.2944064 0.3148678 0.3433022 1.0000000 0.2683965 0.6989141
## Q7_10 0.1445450 0.1264353 0.4485842 0.2334689 0.2683965 1.0000000 0.2664901
## Q7_12 0.5553761 0.3069515 0.2683164 0.3594777 0.6989141 0.2664901 1.0000000
## Q7_14 0.3922082 0.4489550 0.4300696 0.5804649 0.5256589 0.2812001 0.5457326
## Q7_15 0.1813896 0.1093602 0.2558661 0.2474474 0.3159382 0.5161126 0.2604770
## Q7_18 0.2252178 0.3198891 0.2027570 0.3412878 0.2780012 0.2329042 0.2405645
## Q7_21 0.2577641 0.2007712 0.5532761 0.2822111 0.2790424 0.3197975 0.2602531
## Q7_22 0.3994105 0.1852181 0.2325790 0.1978012 0.5694647 0.2151834 0.4634925
## Q7_24 0.2786192 0.3872116 0.3449801 0.5390740 0.2116188 0.2523446 0.2028212
## Q7_14 Q7_15 Q7_18 Q7_21 Q7_22 Q7_24
## Q7_3 0.3922082 0.18138960 0.22521779 0.2577641 0.3994105 0.2786192
## Q7_4 0.4489550 0.10936015 0.31988905 0.2007712 0.1852181 0.3872116
## Q7_5 0.4300696 0.25586613 0.20275702 0.5532761 0.2325790 0.3449801
## Q7_7 0.5804649 0.24744745 0.34128782 0.2822111 0.1978012 0.5390740
## Q7_9 0.5256589 0.31593825 0.27800116 0.2790424 0.5694647 0.2116188
## Q7_10 0.2812001 0.51611265 0.23290420 0.3197975 0.2151834 0.2523446
## Q7_12 0.5457326 0.26047702 0.24056449 0.2602531 0.4634925 0.2028212
## Q7_14 1.0000000 0.21919045 0.42229892 0.3227852 0.2884935 0.4626583
## Q7_15 0.2191904 1.00000000 0.09816021 0.3793122 0.3428206 0.1690617
## Q7_18 0.4222989 0.09816021 1.00000000 0.2490459 0.3142644 0.3727922
## Q7_21 0.3227852 0.37931218 0.24904588 1.0000000 0.3303116 0.2339747
## Q7_22 0.2884935 0.34282056 0.31426445 0.3303116 1.0000000 0.1095123
## Q7_24 0.4626583 0.16906170 0.37279220 0.2339747 0.1095123 1.0000000And high alpha
psych::alpha(s2_bas_final) ##
## Reliability analysis
## Call: psych::alpha(x = s2_bas_final)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.86 0.86 0.89 0.33 6.3 0.017 2.9 0.46 0.31
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.82 0.86 0.89
## Duhachek 0.83 0.86 0.90
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q7_3 0.85 0.85 0.88 0.32 5.7 0.019 0.016 0.30
## Q7_4 0.86 0.86 0.88 0.33 6.0 0.018 0.017 0.30
## Q7_5 0.85 0.85 0.87 0.32 5.7 0.019 0.018 0.29
## Q7_7 0.85 0.85 0.87 0.32 5.6 0.019 0.016 0.29
## Q7_9 0.84 0.85 0.87 0.31 5.5 0.020 0.014 0.30
## Q7_10 0.86 0.86 0.88 0.34 6.1 0.018 0.017 0.32
## Q7_12 0.85 0.85 0.87 0.32 5.6 0.019 0.015 0.31
## Q7_14 0.84 0.84 0.87 0.31 5.4 0.020 0.016 0.28
## Q7_15 0.86 0.86 0.88 0.34 6.2 0.018 0.016 0.32
## Q7_18 0.86 0.86 0.88 0.34 6.1 0.018 0.018 0.32
## Q7_21 0.85 0.86 0.88 0.33 5.9 0.018 0.018 0.31
## Q7_22 0.85 0.86 0.88 0.33 5.9 0.018 0.017 0.31
## Q7_24 0.86 0.86 0.88 0.33 6.0 0.018 0.017 0.31
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q7_3 137 0.64 0.64 0.61 0.56 2.7 0.76
## Q7_4 137 0.54 0.57 0.52 0.47 3.5 0.61
## Q7_5 137 0.65 0.64 0.61 0.55 3.0 0.82
## Q7_7 137 0.67 0.69 0.67 0.59 3.4 0.70
## Q7_9 137 0.73 0.71 0.70 0.65 2.7 0.84
## Q7_10 137 0.55 0.54 0.49 0.45 2.7 0.81
## Q7_12 137 0.69 0.68 0.66 0.61 2.7 0.79
## Q7_14 137 0.73 0.74 0.73 0.67 3.1 0.73
## Q7_15 137 0.52 0.51 0.46 0.42 2.1 0.76
## Q7_18 137 0.53 0.54 0.48 0.43 3.0 0.80
## Q7_21 137 0.59 0.58 0.54 0.50 2.6 0.79
## Q7_22 137 0.60 0.58 0.54 0.50 2.2 0.83
## Q7_24 137 0.54 0.57 0.52 0.47 3.5 0.57
##
## Non missing response frequency for each item
## 1 2 3 4 miss
## Q7_3 0.05 0.31 0.50 0.13 0
## Q7_4 0.00 0.06 0.37 0.57 0
## Q7_5 0.04 0.20 0.46 0.30 0
## Q7_7 0.01 0.10 0.39 0.50 0
## Q7_9 0.08 0.34 0.42 0.16 0
## Q7_10 0.08 0.28 0.50 0.14 0
## Q7_12 0.06 0.31 0.48 0.15 0
## Q7_14 0.02 0.17 0.54 0.27 0
## Q7_15 0.23 0.53 0.21 0.04 0
## Q7_18 0.04 0.21 0.48 0.27 0
## Q7_21 0.07 0.38 0.43 0.12 0
## Q7_22 0.22 0.46 0.26 0.06 0
## Q7_24 0.00 0.04 0.43 0.53 0s2_bas_final$bas_sum <- rowSums(s2_bas_final)
hist(s2_bas_final$bas_sum, breaks=25)
s2_bas_final$ID <- s2_qs$IDModestly inversely related
ComparePars(s2_bas_final$bas_sum, s2_bis_final$bis_sum, use_identity_line = 0)
Prep PTQ s2
Convert to numeric
s2_ptq_num <- foreach (i = 1:ncol(s2_ptq)) %do% {
data.frame(as.numeric(unlist(s2_ptq[i]))) %>% setNames(names(s2_ptq[i]))
}%>% bind_cols()Just 3 data points (< .1% of data) missing for remaining, so mean impute these
length(which(is.na(s2_ptq_num[, 1:15])))## [1] 3length(which(is.na(s2_ptq_num[, 1:15])))/(dim(s2_ptq_num[, 1:15])[1]*dim(s2_ptq_num[, 1:15])[2]) ## [1] 0.001459854s2_ptq_final <- foreach (i = 1:nrow(s2_ptq_num)) %do% {
if (any(is.na(s2_ptq_num[i, ]))) {
s2_ptq_num[i, is.na(s2_ptq_num[i, ])] <- mean(unlist(s2_ptq_num[i, ]), na.rm=TRUE)
}
s2_ptq_num[i, ]
}%>% bind_rows()hist(rowSums(s2_ptq_final[, 1:15]), breaks=25)
s2_ptq_final$ptq_sum <- rowSums(s2_ptq_final[, 1:15])s2_ptq_final$ID <- s2_qs$IDSanity check that correlated in expected direction
ComparePars(s2_ptq_final$ptq_sum, s2_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s2_ptq_final$ptq_sum, s2_bas_final$bas_sum, use_identity_line = 0)
Spot check IDs lined up properly
# s2_ptq_final %>% filter(ID == 22)
# s2_ptq_copy %>% filter(ID == 22)
# # s2_qs %>% filter(ID == 22)
# s2_ptq_final %>% filter(ID == 104)
# s2_ptq_copy %>% filter(ID == 104)
# s2_qs %>% filter(ID == 104)Prep s2 RRS
s2_rrs_num <- foreach (i = 1:ncol(s2_rrs)) %do% {
data.frame(as.numeric(unlist(s2_rrs[i]))) %>% setNames(names(s2_rrs[i]))
}%>% bind_cols()length(which(is.na(s2_rrs_num[, 1:10])))## [1] 3length(which(is.na(s2_rrs_num[, 1:10])))/(dim(s2_rrs_num[, 1:10])[1]*dim(s2_rrs_num[, 1:10])[2]) ## [1] 0.002189781s2_rrs_final <- foreach (i = 1:nrow(s2_rrs_num)) %do% {
if (any(is.na(s2_rrs_num[i, ]))) {
s2_rrs_num[i, is.na(s2_rrs_num[i, ])] <- mean(unlist(s2_rrs_num[i, ]), na.rm=TRUE)
}
s2_rrs_num[i, ]
}%>% bind_rows()hist(rowSums(s2_rrs_final[, 1:10]), breaks=25)
s2_rrs_final$rrs_sum <- rowSums(s2_rrs_final[, 1:10])s2_rrs_final$ID <- s2_qs$IDs2 Spot check
# s2_rrs_final %>% filter(ID == 106)
# s2_rrs_copy %>% filter(ID == 106)
# s2_qs %>% filter(ID == 106) # Q8 is RRSSanity check RRS and PTQ are strongly correlated
ComparePars(s2_rrs_final$rrs_sum, s2_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(s2_rrs_final$rrs_sum, s2_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s2_rrs_final$rrs_sum, s2_bas_final$bas_sum, use_identity_line = 0)
Reduce to just matching IDs
m1_s2_model_red <- m1_study2_eb %>% filter(ID %in% s2_ptq_final$ID)
m1_s2_model_red$c_min_o_phi <- m1_s2_model_red$cog_phi - m1_s2_model_red$overt_phiassert(all(m1_s2_model_red$ID==s2_ptq_final$ID))
assert(all(m1_s2_model_red$ID==s2_rrs_final$ID))PTQ
ComparePars(m1_s2_model_red$c_min_o_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s2_model_red$cog_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)
ComparePars(m1_s2_model_red$overt_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)
Same basic pattern for RRS
ComparePars(m1_s2_model_red$c_min_o_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)
# Would expect higher decay to positively correlate w more brooding — instead small and ns neg
ComparePars(m1_s2_model_red$cog_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)
ComparePars(m1_s2_model_red$overt_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)
s2_rew_train_perf <- s2_train %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_pun_train_perf <- s2_train %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_rew_train_perf$rew_min_pun <- s2_rew_train_perf$m-s2_pun_train_perf$m
s2_rew_sit_perf <- s2_sit %>% filter(valence=="reward") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_pun_sit_perf <- s2_sit %>% filter(valence=="punishment") %>% group_by(ID) %>% summarize(m=mean(correct))
s2_rew_sit_perf$rew_min_pun <- s2_rew_sit_perf$m-s2_pun_sit_perf$m
ComparePars(s2_rew_train_perf$rew_min_pun, s2_rew_sit_perf$rew_min_pun, use_identity_line = 0)
Reduce to just matching IDs
s2_rew_train_perf_red <- s2_rew_train_perf %>% filter(ID %in% s2_bis_final$ID)
s2_rew_sit_perf_red <- s2_rew_sit_perf %>% filter(ID %in% s2_bis_final$ID)
ComparePars(s2_rew_train_perf_red$rew_min_pun, s2_rew_sit_perf_red$rew_min_pun, use_identity_line = 0)
assert(all(s2_rew_sit_perf_red$ID==s2_bis_final$ID))
assert(all(s2_rew_sit_perf_red$ID==s2_bas_final$ID))
assert(all(s2_rew_train_perf_red$ID==s2_bis_final$ID))
assert(all(s2_rew_train_perf_red$ID==s2_bas_final$ID))Not correlated
ComparePars(s2_rew_train_perf_red$rew_min_pun, s2_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s2_rew_train_perf_red$rew_min_pun, s2_bas_final$bas_sum, use_identity_line = 0)
ComparePars(s2_rew_sit_perf_red$rew_min_pun, s2_bis_final$bis_sum, use_identity_line = 0)
ComparePars(s2_rew_sit_perf_red$rew_min_pun, s2_bas_final$bas_sum, use_identity_line = 0)
Added for second revision — Age
There are some Prolific pts without this so take them out and match the datasets
s2_demogs_red <- s2_demogs %>% filter(Age != "DATA_EXPIRED")
s2_demogs_red$Age <- as.numeric(s2_demogs_red$Age)
m1_study2_age_red <- m1_study2_eb %>% filter(ID %in% s2_demogs_red$ID)
assert(all(s2_demogs_red$ID==m1_study2_age_red$ID))m1_study2_age_red$c_min_o_phi <- m1_study2_age_red$cog_phi - m1_study2_age_red$overt_phiOnly one significant is overt decay!
ComparePars(m1_study2_age_red$c_min_o_phi, s2_demogs_red$Age, use_identity_line = 0)
ComparePars(m1_study2_age_red$cog_phi, s2_demogs_red$Age, use_identity_line = 0)
ComparePars(m1_study2_age_red$overt_phi, s2_demogs_red$Age, use_identity_line = 0)
Age and task performance
s2_train_perf_age_red <- s2_train %>% filter(ID %in% s2_demogs_red$ID) %>% group_by(ID) %>% summarize(m=mean(correct))
s2_test_perf_age_red <- s2_sit %>% filter(ID %in% s2_demogs_red$ID) %>% group_by(ID) %>% summarize(m=mean(correct))assert(all(s2_demogs_red$ID==s2_train_perf_age_red$ID))
assert(all(s2_demogs_red$ID==s2_test_perf_age_red$ID))In study 2 unrelated to overall performance
ComparePars(s2_train_perf_age_red$m, s2_demogs_red$Age, use_identity_line = 0)
ComparePars(s2_test_perf_age_red$m, s2_demogs_red$Age, use_identity_line = 0)
s2_perf_age_red <- s2_train %>% filter(ID %in% s2_demogs_red$ID) %>%
group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_perf_s2 <- as.numeric(unlist(s2_perf_age_red[s2_perf_age_red$condition=="overt","m"]-
s2_perf_age_red[s2_perf_age_red$condition=="cognitive", "m"]))
# .. and test subject-level differences in performance
s2_perf_age_red_test <- s2_sit %>% filter(ID %in% s2_demogs_red$ID) %>%
group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_testperf_s2 <- as.numeric(unlist(s2_perf_age_red_test[s2_perf_age_red_test$condition=="overt","m"]-
s2_perf_age_red_test[s2_perf_age_red_test$condition=="cognitive", "m"]))ComparePars(o_min_c_perf_s2, s2_demogs_red$Age, use_identity_line = 0)
ComparePars(o_min_c_testperf_s2, s2_demogs_red$Age, use_identity_line = 0)